Iterative solver algorithm
Iterative solvers algorithms are described in a very extensive way in the textbook of Yousef Saad [1]. Probably the first iterative solver was the method of solving the 4 by 4 equation system proposed by Gauss. Discussing the general family of iterative solver algorithms is beyond the scope of this manual. In this module, we will discuss an iterative solver applicable to the projection problem where the area has any shape that is not a regular square.
Our example problem to be solved is a two-dimensional projection problem defined in an area that does not have a square shape (for example, the problem of projecting a bitmap onto a folded area). Then our system of equations to solve looks like this \( \begin{bmatrix} \int_{\Omega}{B^x_{1,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int_{\Omega}{B^x_{1,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{\Omega}{B^x_{1,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \int_{\Omega}{B^x_{2,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int_{\Omega}{B^x_{2,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{\Omega}{B^x_{2,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \vdots & \vdots & \vdots & \vdots \\ \int_{\Omega}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{1,p}B^y_{1,p}} & \int_{\Omega}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{\Omega}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \end{bmatrix} \begin{bmatrix} {\color{red}u_{1,1}} \\ {\color{red}u_{2,1}} \\ \vdots \\ {\color{red}u_{N_x,N_y}}\\ \end{bmatrix} \)
\( = \begin{bmatrix} \int_{\Omega} BITMAP(x,y) {\color{blue}B^x_1(x)*B^y_1(y)} \\ \int_{\Omega} BITMAP(x,y) {\color{blue}B^x_1(x)*B^y_2(y)} \\ \vdots \\ \int_{\Omega} BITMAP(x,y) {\color{blue}B^x_{N_x}(x)*B^y_{N_y}(y)} \\ \end{bmatrix} \)
where \( \Omega \) where Ω is no longer a square but our "folded" area.
The derivation of this system of equations proceeds in the same way as described in chapter one for bitmap projection tasks, the only difference is that we are now projecting onto a "folded" area. For example, we want to calculate and visualize what a bitmap will look like when superimposed over a pleated material (e.g., a curtain).
Of course, for each integral over the area of \( \Omega \) we can change the variables to a square area, and then the jacobian of this transformation will be included in the integral
\( \int_{\Omega} {\color{blue}B^x_{i}(x)B^y_{j}(y)} \color{red}{B^x_{k}(x)B^y_{l}(y)} dxdy = \int_{[0,1]^2} \color{blue}{B^x_{i}(x)B^y_{j}(y)} \color{red}{B^x_{k}(x)B^y_{l}(y)} |Jac(x,y)| dxdy \)
We would like to use a variable-direction solver to factorize our matrix, however, this is not possible because we cannot separate variables in Jacobian, we do not know how to process \( |Jac(x,y)|=|Jac_x(x)|*|Jac_y(y)| \).
Fortunately, we can apply the following preconditioner iterative algorithm
Our goal is to solve a system of equations \( Ax=y \), where \( A \) is, for example, the matrix of our bitmap projection problem into a "rippled" area. The size of the matrix \( A \) to \( N_xN_y \) while its terms are given in our example by \( {A}_{i,j}=\int_{[0,1]^2} \color{blue}{B^x_{i}(x)B^y_{j}(y)} \color{red}{B^x_{k}(x)B^y_{l}(y)} |Jac(x,y)| dx \).
- We define our "preconditioner" or matrix \( M \) which arises from the matrix \( A \) by settig the Jacobian to one. In other words, the matrix terms are given by \( {M}_{i,j}=\int_{[0,1]^2}\color{blue}{B^x_{i}(x)B^y_{j}(y)} \color{red}{B^x_{k}(x)B^y_{l}(y)} dx \). In our example, this means replacing the projection into an undulating area with a projection into a flat area.
- Let as denote \( y=Mx \)
- We choose any "starting point" (any vector \( y_0 \) B-spline approximation coefficients of \( N_xN_y \)). We can take for example \( {y_0}_{i,j}= \) pixel value at the point where the maximum of the B-spline value is \( B^x_i(x) B^y_j(y) \). But I might as well take \( {y_0} _ {i, j} = 0 \) and our algorithm will also work properly (although it will probably require a bit more iteration).
- Then we calculate the so-called residuum \( r= M^{-1} b - M^{-1} A M^{-1} y_0 = M^{-1} (b - A M^{-1} M y_0) \). We do it in the following steps
- denote \( M^{-1}y_0=z_0 \) and we solve a system of equations \( Mz_0=y_0 \). Note that this means solving the same system of equations as in chapter one.
\( \begin{bmatrix} \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \vdots & \vdots & \vdots & \vdots \\ \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \end{bmatrix} \begin{bmatrix} {\color{red}{z_0}_{1,1}} \\ {\color{red}{z_0}_{2,1}} \\ \vdots \\ {\color{red}{z_0}_{N_x,N_y}}\\ \end{bmatrix}= \) \( \begin{bmatrix} {y_0}_{1,1} \\ {y_0}_{2,1} \\ \vdots \\ {y_0}_{N_x,N_y} \end{bmatrix} \)
We can do this with a linear-directional variable complexity.
- Then we have \( r=M^{-1}(b-Az_0) \) where \( b-Az_0=v \) this is the solution coefficient vector, i.e. \( M^{-1}v=r \). We need to multiply \( A \) (with Jacobian) by the resulting solution coefficient vector. We're going to get a vector
\( v=b-Az_0 = \begin{bmatrix} v_{1,1} \\ v_{2,1} \\ \cdots \\ v_{N_x,N_y} \\ \end{bmatrix} = \\ \begin{bmatrix} \int_{\Omega} BITMAP(x,y) {\color{blue}B^x_1(x)*B^y_1(y)} - \sum_{i=1,...,N_x; j=1,...,N_y }{z_0}_{i,j} \int_{\Omega} B^x_{1,p} B^y_{1,p} B^x_{i,p}B^y_{j,p} \\ \int_{\Omega} BITMAP(x,y) {\color{blue}B^x_2(x)*B^y_1(y) } - \sum_{i=1,...,N_x; j=1,...,N_y} {z_0}_{i,j } \int_{\Omega} B^x_{2,p }B^y_{1,p} B^x_{i,p}B^y_{j,p} \\ \cdots \\ \int_{\Omega} BITMAP(x,y) {\color{blue}B^x_{N_x}(x)*B^y_{N_y}(y) } - \sum_{i=1,...,N_x; j=1,...,N_y }{z_0}_{i,j} \int_{\Omega } B^x_{N_x,p}B^y_{N_y,p} B^x_{i,p}B^y_{j,p} \\ \end{bmatrix} \)
- We solve \( Mr=v \).
\( \begin{bmatrix} \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \vdots & \vdots & \vdots & \vdots \\ \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \end{bmatrix} \)
\( \begin{bmatrix} {\color{red}r_{1,1}} \\ {\color{red}r_{2,1}} \\ \vdots \\ {\color{red}r_{N_x,N_y}}\\ \end{bmatrix} \)
\( = \begin{bmatrix} \int_{\Omega} BITMAP(x,y) {\color{blue}B^x_1(x)*B^y_1(y)} - \sum_{i=1,...,N_x;j=1,...,N_y}{z_0}_{i,j}\int_{\Omega}{B^x_{1,p}B^y_{1,p}}{B^x_{i,p}B^y_{j,p}} \\ \int_{\Omega } BITMAP(x,y) {\color{blue}B^x_2(x)*B^y_1(y)} - \sum_{i=1,...,N_x;j=1,...,N_y}{z_0}_{i,j}\int_{\Omega}B^x_{2,p}B^y_{1,p}B^x_{i,p}B^y_{j,p} \\ \cdots \\ \int_{\Omega } BITMAP(x,y) {\color{blue}B^x_{N_x}(x)*B^y_{N_y}(y)} - \sum_{i=1,...,N_x;j=1,...,N_y}{z_0}_{i,j}\int_{\Omega}B^x_{N_x,p}B^y_{N_y,p}B^x_{i,p}B^y_{j,p } \\ \end{bmatrix} \)
Again, we can do this with a alternating-direction solver with linear computational complexity.
Having the residuum counted \( r \) we can correct our proposed solution \( y_1=y_0+r \)
The melt condition is sufficiently low residuum. In other words, I have to count the quota from the residuum (check how small the whole residuum is, for example \( ||r||=\left( \sum_{i=1,...,Nx;j=1,...,Ny} |r_{i,j}|^2\right)^{0.5} \) If \( ||r||>\epsilon \) we mark \( y_0:=y_1 \) and go to the point 5.
Finally, we calculate our final solution (bitmap projection coefficients on the corrugated area) by solving the system of equations \( Mx_{final}=y_{final} \). Our final solution is \( y_{final} \).
\( \begin{bmatrix} \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \vdots & \vdots & \vdots & \vdots \\ \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \end{bmatrix} \begin{bmatrix} {\color{red}{x_{final}}_{1,1}} \\ {\color{red}{x_{final}}_{2,1}} \\ \vdots \\ {\color{red}{x_{final}}_{N_x,N_y}}\\ \end{bmatrix} =\\ = \begin{bmatrix} {y_{final}}_{1,1} \\ {y_{final}}_{2,1} \\ \vdots \\ {y_{final}}_{N_x,N_y} \\ \end{bmatrix} \)
Again, we can do this with an alternating-direction solver with linear computational complexity.